rm(list=ls())

root = "/Volumes/wcnr-network/Research/Hobbs/A_Bison_2011/Analysis for Ecological Monograph/"
root2="/Volumes/wcnr-network/Research/Hobbs/A_Bison_2011/Analysis for Ecological Monograph/Figures and Tables/Results/"

#root = "/Users/Tom/Documents/YNP_Bison_Model/2011 Model Update/Analysis for Manuscript/"
#root2="/Users/Tom/Documents/YNP_Bison_Model/A_Bison_Manuscript/"

path = function (stem,r=root){
return(paste(r,stem, sep=""))	
}

library(rjags)
library(coda)
library(xtable)

setwd(path("Model and analysis code/Analysis/"))

load(path("Data and data reduction/Data objects/2011 Bison Data with multinomial.Rdata"))

pdf(file=path("All_data.pdf", r=root2), height=7,width=4)


y.errorbars = function(x,y,ybar, ylab=x, xlab=y, main=" ",y.low,y.high){
	plot(x,y,pch=16, cex=.6, xlab=xlab, ylab = ylab, main=main, ylim=c(y.low,y.high), family="serif")
	arrows(x,y-ybar,x,y+ybar, code=3, angle=90, length=.01, lwd=.5)
}

par(mfrow=c(4,2))
#====census data
w=data$census$N
par(cex=.6)
y.errorbars(x=w$year,y=w$count.mean,ybar=2*w$count.sd, y.low=0, y.high=6000, xlab="Year", ylab="Population count", main="A")
points(data$removals$count$year,data$removals$count$removal,typ="h")

#=====aerial calf data
w=data$aerial$calf
y.errorbars(x=w$year,y=w$mean,ybar=2*w$sd_of_mean, y.low=0, y.high=.6, xlab="Year", ylab="Proportion of population", main="B")

#=====ground yearlings
w=as.data.frame(data$ground)
y.errorbars(x=w$mu.ground.year,y=w$mu.yrlg_mean,ybar=2*w$sigma.yrlg_sd, y.low=0, y.high=.6, xlab="Year", ylab="Proportion of population", main="C")
#=====ground cowsw=as.data.frame(data$ground)
y.errorbars(x=w$mu.ground.year,y=w$mu.cow_mean,ybar=2*w$sigma.cow_sd, y.low=0, y.high=.6, xlab="Year", ylab="Proportion of population", main="D")

#=====ground cows
w=as.data.frame(data$ground)
y.errorbars(x=w$mu.ground.year,y=w$mu.bull_mean,ybar=2*w$sigma.bull_sd, y.low=0, y.high=.6, xlab="Year", ylab="Proportion of population", main="E")

####calf serology
a=data$sero$calf[,3]+1
b=data$sero$calf[,4]-data$sero$calf[,3]+1
mu=a/(a+b)
stdev = sqrt(a*b/((a=b^2)*(a+b+1)))
y.data=as.data.frame(cbind(data$sero$calf,mu,stdev))
year=as.data.frame(1989:2010)
names(year)=c("year")
w=merge(y.data,year,by="year",all.y=TRUE)


y.errorbars(x=w$year,y=w$mu,ybar=2*w$stdev, y.low=0, y.high=1, xlab="Year", ylab="Proportion seropositive", main="F")

####yrlg serology
a=data$sero$yrl[,3]+1
b=data$sero$yrl[,4]-data$sero$yrl[,3]+1
mu=a/(a+b)
stdev = sqrt(a*b/((a=b^2)*(a+b+1)))
y.data=as.data.frame(cbind(data$sero$yrl,mu,stdev))
year=as.data.frame(1989:2010)
names(year)=c("year")
w=merge(y.data,year,by="year",all.y=TRUE)


y.errorbars(x=w$year,y=w$mu,ybar=2*w$stdev, y.low=0, y.high=1, xlab="Year", ylab="Proportion seropositive", main="G")

####cow serology
a=data$sero$cow[,3]+1
b=data$sero$cow[,4]-data$sero$cow[,3]+1
mu=a/(a+b)
stdev = sqrt(a*b/((a=b^2)*(a+b+1)))
y.data=as.data.frame(cbind(data$sero$cow,mu,stdev))
year=as.data.frame(1989:2010)
names(year)=c("year")
w=merge(y.data,year,by="year",all.y=TRUE)


y.errorbars(x=w$year,y=w$mu,ybar=2*w$stdev, y.low=0, y.high=1, xlab="Year", ylab="Proportion seropositive", main="H")

dev.off()